#################################
##' Code for:
##' Title: "Effective climate clubs require ambition, leverage, and insulation"
##' Author: Sam S. Rowan
##' Contact: sam.rowan [at] concordia.ca
##' Date: 2023-10-19
#################################
##' This script:
##' Dropping the most constrained core member when setting the trade penalty
#################################

#### Initialize environment ####

## Packages
library(tidyverse)
library(truncnorm)
library(knitr)

## Working directory
# setwd("/Users/srowa/Dropbox/Projects/climate_clubs_trade")
# -- Set working directory to replication folder

## Define the logistic function
logistic_function <- function(x, k, a) {
  return(1 / (1 + exp(-k * (x - a))))
}

## Plot settings for later
theme_clubs <-  theme(legend.position = "bottom", 
                      panel.background = element_rect(fill=NA),
                      panel.grid.major.x = element_blank(),
                      panel.grid.major.y = element_blank(),
                      panel.grid.minor.y =  element_blank(),
                      panel.grid.minor.x =  element_blank(),
                      panel.ontop = FALSE,
                      axis.line = element_line(linewidth = 1, colour = "black"),
                      axis.text = element_text(size = 16),
                      legend.text = element_text(size = 12),
                      legend.title = element_text(size = 16),
                      axis.title = element_text(size = 16))


#### Define climate model parameters ####

## Define number of states
num_states <- 150

## Draw their climate ideal points
set.seed(202109)
climate_ideal <- cbind.data.frame(
  climate_ideal = runif(n = num_states, min = 0, max = 1)) %>% 
  arrange(climate_ideal) %>% 
  mutate(id = seq(from = 1, to = num_states, by = 1))
# climate_ideal

####' ---- NOTE: Code block below runs for about 5'
####' ---- The simulation output is saved as 'output/members_drop_constrained_theta75.csv'
####' ---- Skip running the simulation and load this in directly at around line 275
####' ---- And build figures from there


#### Now simulate this 500 times ####

## Create loop counters
set.seed(202109)
seeds_vector <- as.integer(runif(n = 500, min = 1, max = 1e6))
theta_vector <- 0.75

## Create objects to hold outputs
# Count membership for each value of theta and each draw of trade ties
density_members_independent <- matrix(NA, nrow = length(theta_vector), ncol = length(seeds_vector))
density_members_positive <- matrix(NA, nrow = length(theta_vector), ncol = length(seeds_vector))
density_members_negative <- matrix(NA, nrow = length(theta_vector), ncol = length(seeds_vector))

# Takes about 5'
start_time <- Sys.time()
for (seed in seeds_vector) {
  # Declare seed
  set.seed(seed)
  theta_df <-  cbind.data.frame(
    # Add IDs
    id_a = rep(climate_ideal$id, times = num_states),
    id_b = rep(climate_ideal$id, each = num_states),
    # Bind climate ideal points into dyadic setup
    climate_ideal_a = rep(climate_ideal$climate_ideal, times = num_states),
    climate_ideal_b = rep(climate_ideal$climate_ideal, each = num_states)) %>%
    # Measure dyadic ideal point distance
    mutate(ideal_point_distance = abs(climate_ideal_a - climate_ideal_b)) %>%
    # Draw trade data based on ideal point distance
    mutate(trade_independent = rtruncnorm(n = num_states^2,
                                          # a = 0, b = 1,
                                          mean = 0,
                                          sd = 0.5),
           trade_positive = rtruncnorm(n = num_states^2,
                                       # a = 0, b = 1,
                                       mean = (1-ideal_point_distance), # negative of difference, implies trade more when ideal points closer
                                       sd = 0.5),
           trade_negative = rtruncnorm(n = num_states^2,
                                       # a = 0, b = 1,
                                       mean = ideal_point_distance, # positive of difference, implies trade less when ideal points closer
                                       sd = 0.5)) %>% 
    # Use logistic function to bound trade ties
    mutate_at(vars(trade_independent:trade_negative),
              ~logistic_function(x = ., k = 10, a = 0.75)) %>% 
    ## Calculate dyadic trade flows as shares
    # Remove data for faux dyads -- where id_a == id_b
    mutate_at(vars(ideal_point_distance,
                   trade_independent, trade_positive, trade_negative),
              ~replace(., id_a == id_b, NA)) %>% 
    # Get each state's total trade
    group_by(id_a) %>% 
    mutate(sum_trade_independent = sum(trade_independent, na.rm = T),
           sum_trade_positive = sum(trade_positive, na.rm = T),
           sum_trade_negative = sum(trade_negative, na.rm = T)) %>% 
    # Standardize
    rowwise() %>% 
    mutate(share_dyadic_independent = trade_independent / sum_trade_independent,
           share_dyadic_positive = trade_positive / sum_trade_positive,
           share_dyadic_negative = trade_negative / sum_trade_negative) %>% 
    select(id_a, id_b, climate_ideal_a, climate_ideal_b, ideal_point_distance,
           share_dyadic_independent, share_dyadic_positive, share_dyadic_negative) %>% 
    ## Calculate each state's utility to membership
    # Create a club's common climate goal
    mutate(theta = theta_vector,
           core_member_a = ifelse(climate_ideal_a >= theta, 1, 0),
           core_member_b = ifelse(climate_ideal_b >= theta, 1, 0)) %>% 
    # Flag trade to club members
    mutate(trade_to_club_independent = ifelse(core_member_b == 1, 
                                              share_dyadic_independent, NA),
           trade_to_club_positive = ifelse(core_member_b == 1, 
                                           share_dyadic_positive, NA),
           trade_to_club_negative = ifelse(core_member_b == 1, 
                                           share_dyadic_negative, NA),) %>% 
    # Get share of trade to all club members
    group_by(id_a) %>% 
    mutate(share_trade_club_independent = sum(trade_to_club_independent, na.rm = T),
           share_trade_club_positive = sum(trade_to_club_positive, na.rm = T),
           share_trade_club_negative = sum(trade_to_club_negative, na.rm = T)) %>% 
    ungroup() %>% 
    # Identify most constrained and second most constrained members
    group_by(core_member_a) %>% 
    mutate(most_constrained_independent = dplyr::nth(x = share_trade_club_independent,
                                                   n = 1L, 
                                                   order_by = share_trade_club_independent),
         most_constrained_positive = dplyr::nth(x = share_trade_club_positive,
                                                n = 1L, 
                                                order_by = share_trade_club_positive),
         most_constrained_negative = dplyr::nth(x = share_trade_club_negative,
                                                n = 1L, 
                                                order_by = share_trade_club_negative),
         # The 151th element is after all 150 dyadic flows for the most constrained state
         second_constrained_independent = dplyr::nth(x = share_trade_club_independent,
                                                     n = 151L, 
                                                     order_by = share_trade_club_independent),
         second_constrained_positive = dplyr::nth(x = share_trade_club_positive,
                                                  n = 151L, 
                                                  order_by = share_trade_club_positive),
         second_constrained_negative = dplyr::nth(x = share_trade_club_negative,
                                                  n = 151L, 
                                                  order_by = share_trade_club_negative)) %>%
    ungroup() %>% 
    mutate_at(vars(most_constrained_independent:second_constrained_negative),
              ~ifelse(core_member_a == 0, NA, .)) %>% 
    # Overwrite their club membership
    mutate(core_member_a_independent = ifelse(
      share_trade_club_independent == 
        most_constrained_independent, 0, core_member_a),
      core_member_b_independent = ifelse(
        share_trade_club_independent == 
          most_constrained_independent, 0, core_member_b),
      core_member_a_positive = ifelse(
        share_trade_club_positive == 
          most_constrained_positive, 0, core_member_a),
      core_member_b_positive = ifelse(
        share_trade_club_positive == 
          most_constrained_positive, 0, core_member_b),
      core_member_a_negative = ifelse(
        share_trade_club_negative == 
          most_constrained_negative, 0, core_member_a),
      core_member_b_negative = ifelse(
        share_trade_club_negative == 
          most_constrained_negative, 0, core_member_b)) %>% 
    mutate_at(vars(core_member_a_independent:core_member_b_negative),
              ~ifelse(is.na(.), 0, .)) %>% 
    # Now re-calculate the trade to club variables 
    mutate(trade_two_club_independent = ifelse(core_member_b == 1, 
                                               share_dyadic_independent, NA),
           trade_two_club_positive = ifelse(core_member_b == 1, 
                                            share_dyadic_positive, NA),
           trade_two_club_negative = ifelse(core_member_b == 1, 
                                            share_dyadic_negative, NA)) %>% 
    # Get share of trade to all club members
    group_by(id_a) %>% 
    mutate(share_trade_two_club_independent = sum(trade_two_club_independent, na.rm = T),
           share_trade_two_club_positive = sum(trade_two_club_positive, na.rm = T),
           share_trade_two_club_negative = sum(trade_two_club_negative, na.rm = T)) %>% 
    ungroup() %>% 
    # Define lambda as second most constrained
    mutate(lambda_two_independent = min(second_constrained_independent, na.rm = T),
           lambda_two_positive = min(second_constrained_positive, na.rm = T),
           lambda_two_negative = min(second_constrained_negative, na.rm = T)) %>% 
    # Calculate change in lambda 
    mutate(delta_lambda_independent = second_constrained_independent - most_constrained_independent,
           delta_lambda_positive = second_constrained_positive - most_constrained_positive,
           delta_lambda_negative = second_constrained_negative - most_constrained_negative) %>% 
    distinct(climate_ideal_a, theta, share_trade_two_club_independent, 
             share_trade_two_club_positive, share_trade_two_club_negative,
             lambda_two_independent, lambda_two_positive, lambda_two_negative,
             delta_lambda_independent, delta_lambda_positive, delta_lambda_negative, ) %>% 
    # Calculate utilities to membership
    mutate(beta = 1) %>% # Define weights for climate versus trade
    rowwise() %>% 
    # Calculate utilities for the *independent* world
    # Utility to join
    mutate(utility_join_independent = (-abs(climate_ideal_a - theta)) + # Climate
             (1 - share_trade_two_club_independent) * (1 - lambda_two_independent)*beta + # Non-club trade
             (share_trade_two_club_independent*beta) - # Club trade
             (theta - lambda_two_independent)) %>%  # Aversion
    # Utility to abstain
    mutate(utility_abstain_independent = (-climate_ideal_a) + # Climate
             (1-share_trade_two_club_independent)*beta + # Non-club trade
             (share_trade_two_club_independent*(1-lambda_two_independent)*beta) + # Club trade
             (0)) %>%  # Aversion
    # Net, membership
    mutate(net_utility_independent = utility_join_independent - utility_abstain_independent,
           membership_independent = ifelse(net_utility_independent >= 0, 1, 0)) %>% 
    # Calculate utilities for the *positive* world
    # Utility to join
    mutate(utility_join_positive = (-abs(climate_ideal_a - theta)) + # Climate
             (1 - share_trade_two_club_positive) * (1 - lambda_two_positive)*beta + # Non-club trade
             (share_trade_two_club_positive*beta) - # Club trade
             (theta - lambda_two_positive)) %>%  # Aversion
    # Utility to abstain
    mutate(utility_abstain_positive = (-climate_ideal_a) + # Climate
             (1-share_trade_two_club_positive)*beta + # Non-club trade
             (share_trade_two_club_positive*(1-lambda_two_positive)*beta) + # Club trade
             (0)) %>%  # Aversion
    # Net, membership
    mutate(net_utility_positive = utility_join_positive - utility_abstain_positive,
           membership_positive = ifelse(net_utility_positive >= 0, 1, 0)) %>% 
    # Calculate utilities for the *negative* world
    # Utility to join
    mutate(utility_join_negative = (-abs(climate_ideal_a - theta)) + # Climate
             (1 - share_trade_two_club_negative) * (1 - lambda_two_negative)*beta + # Non-club trade
             (share_trade_two_club_negative*beta) - # Club trade
             (theta - lambda_two_negative)) %>%  # Aversion
    # Utility to abstain
    mutate(utility_abstain_negative = (-climate_ideal_a) + # Climate
             (1-share_trade_two_club_negative)*beta + # Non-club trade
             (share_trade_two_club_negative*(1-lambda_two_negative)*beta) + # Club trade
             (0)) %>%  # Aversion
    # Net, membership
    mutate(net_utility_negative = utility_join_negative - utility_abstain_negative,
           membership_negative = ifelse(net_utility_negative >= 0, 1, 0))
  
  # Save the number of countries in the club
  density_members_independent[which(seed == seeds_vector)] <- nrow(
    filter(theta_df, membership_independent == 1))
  density_members_positive[which(seed == seeds_vector)] <- nrow(
    filter(theta_df, membership_positive == 1))
  density_members_negative[which(seed == seeds_vector)] <- nrow(
    filter(theta_df, membership_negative == 1))
}
end_time <- Sys.time()
end_time - start_time

## Summarize outputs
members_constrained <- cbind.data.frame(independent = density_members_independent[1,],
                 positive = density_members_positive[1,],
                 negative = density_members_negative[1,]) %>% 
  as.data.frame() %>% 
  mutate(run = "Drop most constrained core member")

## Write to disk
# write_csv(x = members_constrained,
#           file = "output/members_drop_constrained_theta75.csv")
members_constrained <- read_csv("output/members_drop_constrained_theta75.csv")

## Not run --- some summaries
# # Table
# members_constrained %>% 
#   pivot_longer(names_to = "model", 
#                values_to = "members",
#                independent:negative) %>% 
#   group_by(model) %>% 
#   summarize(median = median(members)) %>% 
#   mutate(prop.median = median/150)
# 
# # Plot
# members_constrained %>% 
#   pivot_longer(names_to = "model", 
#                values_to = "members",
#                independent:negative) %>% 
#   group_by(model, members) %>% 
#   summarize(n = n()) %>% 
#   ggplot(., aes(x = members, y = n, fill = model)) +
#   geom_col(position = "dodge") +
#   scale_fill_manual(values = c("#7DCEA0", "#229954", "#556B2F")) +
#   labs(x = "Members", y = "Number of simulations", fill = "Model") +
#   theme_clubs 


#### Compare to baseline model ####

## Load membership summaries from baseline simulation

members_baseline <- read_csv("output/members_baseline_theta75.csv")

## Summary of how many additional members
bind_rows(members_baseline,
          members_constrained) %>% 
  mutate(model = ifelse(str_detect(run, "constrained"), "constrained", "baseline")) %>% 
  mutate_at(vars(independent:negative), ~as.numeric(.)) %>% 
  mutate(run = rep(seq(1:500), times = 2)) %>% 
  pivot_longer(names_to = "correlation",
               values_to = "members",
               independent:negative) %>% 
  pivot_wider(names_from = model,
              values_from = members) %>% 
  mutate(difference = constrained - baseline) %>% 
  group_by(correlation) %>% 
  summarize(mean_difference = mean(difference))

## Combine outputs and plot
additional_members_at_theta75 <- bind_rows(members_baseline,
          members_constrained) %>% 
  mutate(model = ifelse(str_detect(run, "constrained"), "constrained", "baseline")) %>% 
  mutate_at(vars(independent:negative), ~as.numeric(.)) %>% 
  mutate(run = rep(seq(1:500), times = 2)) %>% 
  pivot_longer(names_to = "correlation",
               values_to = "members",
               independent:negative) %>% 
  pivot_wider(names_from = model,
              values_from = members) %>% 
  mutate(difference = constrained - baseline) %>% 
  select(run, correlation, difference) %>% 
  group_by(correlation, difference) %>% 
  summarize(n = n()) %>% 
  mutate(model = case_when(correlation == "independent" ~ "Independent",
                           correlation == "positive" ~ "Positive",
                           correlation == "negative" ~ "Negative")) %>% 
  ggplot(., aes(x = difference, y = n, fill = model)) +
  # geom_col(position = "dodge") +
  geom_col(position = position_dodge2(width = 0.9, preserve = "single")) +
  scale_x_continuous(breaks = seq(0, 10, 1)) +
  scale_fill_manual(values = c("#7DCEA0", "#229954", "#556B2F")) +
  labs(x = expression(Additional~members~at~theta[j]==0.75),
       y = "Number of simulations",
       fill = "Model") +
  theme_clubs
additional_members_at_theta75
ggsave(plot = additional_members_at_theta75, filename = "text/additional_members_at_theta75.pdf", units = 'in', height = 6, width = 8)




#### NOT RUN: One example to show how lambda changes ####

#' ## Draw dyadic trade flows as a function of ideal points
#' climate_trade_data <- cbind.data.frame(
#'   # Add IDs
#'   id_a = rep(climate_ideal$id, times = num_states),
#'   id_b = rep(climate_ideal$id, each = num_states),
#'   # Bind climate ideal points into dyadic setup
#'   climate_ideal_a = rep(climate_ideal$climate_ideal, times = num_states),
#'   climate_ideal_b = rep(climate_ideal$climate_ideal, each = num_states)) %>%
#'   # Measure dyadic ideal point distance
#'   mutate(ideal_point_distance = abs(climate_ideal_a - climate_ideal_b)) %>%
#'   # Draw trade data based on ideal point distance
#'   mutate(trade_independent = rtruncnorm(n = num_states^2,
#'                                         # a = 0, b = 1,
#'                                         mean = 0,
#'                                         sd = 0.5),
#'          trade_positive = rtruncnorm(n = num_states^2,
#'                                      # a = 0, b = 1,
#'                                      mean = (1-ideal_point_distance), # negative of difference, implies trade more when ideal points closer
#'                                      sd = 0.5),
#'          trade_negative = rtruncnorm(n = num_states^2,
#'                                      # a = 0, b = 1,
#'                                      mean = ideal_point_distance, # positive of difference, implies trade less when ideal points closer
#'                                      sd = 0.5)) %>% 
#'   # Use logistic function to bound trade ties
#'   mutate_at(vars(trade_independent:trade_negative),
#'             ~logistic_function(x = ., k = 10, a = 0.75))
#' 
#' ## Calculate dyadic trade flows as shares
#' climate_trade_shares <- climate_trade_data %>% 
#'   # Remove data for faux dyads -- where id_a == id_b
#'   mutate_at(vars(ideal_point_distance,
#'                  trade_independent, trade_positive, trade_negative),
#'             ~replace(., id_a == id_b, NA)) %>% 
#'   # Get each state's total trade
#'   group_by(id_a) %>% 
#'   mutate(sum_trade_independent = sum(trade_independent, na.rm = T),
#'          sum_trade_positive = sum(trade_positive, na.rm = T),
#'          sum_trade_negative = sum(trade_negative, na.rm = T)) %>% 
#'   # Standardize
#'   rowwise() %>% 
#'   mutate(share_dyadic_independent = trade_independent / sum_trade_independent,
#'          share_dyadic_positive = trade_positive / sum_trade_positive,
#'          share_dyadic_negative = trade_negative / sum_trade_negative) %>% 
#'   select(id_a, id_b, climate_ideal_a, climate_ideal_b, ideal_point_distance,
#'          share_dyadic_independent, share_dyadic_positive, share_dyadic_negative)
#' #' -- Interpretation
#' #' Each state has their own climate ideal point _a _b
#' #' The trade variables are:
#' #' How much _a exports to _b
#' #' share_dyadic_independent: what share of _a's exports are to _b
#' 
#' ## Define a value of theta
#' theta_vector <- 0.75
#' 
#' ## Calculate each state's utility to membership
#' climate_trade_utilities <- climate_trade_shares %>%
#'   # Create a club's common climate goal
#'   mutate(theta = theta_vector,
#'          core_member_a = ifelse(climate_ideal_a >= theta, 1, 0),
#'          core_member_b = ifelse(climate_ideal_b >= theta, 1, 0)) %>% 
#'   # Flag trade to club members
#'   mutate(trade_to_club_independent = ifelse(core_member_b == 1, 
#'                                             share_dyadic_independent, NA),
#'          trade_to_club_positive = ifelse(core_member_b == 1, 
#'                                          share_dyadic_positive, NA),
#'          trade_to_club_negative = ifelse(core_member_b == 1, 
#'                                          share_dyadic_negative, NA),) %>% 
#'   # Get share of trade to all club members
#'   group_by(id_a) %>% 
#'   mutate(share_trade_club_independent = sum(trade_to_club_independent, na.rm = T),
#'          share_trade_club_positive = sum(trade_to_club_positive, na.rm = T),
#'          share_trade_club_negative = sum(trade_to_club_negative, na.rm = T)) %>% 
#'   ungroup() %>% 
#'   # Identify most constrained and second most constrained members
#'   group_by(core_member_a) %>% 
#    mutate(most_constrained_independent = dplyr::nth(x = share_trade_club_independent,
#                                                  n = 1L, 
#                                                  order_by = share_trade_club_independent),
#        most_constrained_positive = dplyr::nth(x = share_trade_club_positive,
#                                               n = 1L, 
#                                               order_by = share_trade_club_positive),
#        most_constrained_negative = dplyr::nth(x = share_trade_club_negative,
#                                               n = 1L, 
#                                               order_by = share_trade_club_negative),
#        # The 151th element is after all 150 dyadic flows for most constrained
#        second_constrained_independent = dplyr::nth(x = share_trade_club_independent,
#                                                    n = 151L, 
#                                                    order_by = share_trade_club_independent),
#        second_constrained_positive = dplyr::nth(x = share_trade_club_positive,
#                                                 n = 151L, 
#                                                 order_by = share_trade_club_positive),
#        second_constrained_negative = dplyr::nth(x = share_trade_club_negative,
#                                                 n = 151L, 
#                                                 order_by = share_trade_club_negative)) %>%
#'   ungroup() %>% 
#'   mutate_at(vars(most_constrained_independent:second_constrained_negative),
#'             ~ifelse(core_member_a == 0, NA, .)) %>% 
#'   # Overwrite their club membership
#'   mutate(core_member_a_independent = ifelse(
#'     share_trade_club_independent == 
#'       most_constrained_independent, 0, core_member_a),
#'     core_member_b_independent = ifelse(
#'       share_trade_club_independent == 
#'         most_constrained_independent, 0, core_member_b),
#'     core_member_a_positive = ifelse(
#'       share_trade_club_positive == 
#'         most_constrained_positive, 0, core_member_a),
#'     core_member_b_positive = ifelse(
#'       share_trade_club_positive == 
#'         most_constrained_positive, 0, core_member_b),
#'     core_member_a_negative = ifelse(
#'       share_trade_club_negative == 
#'         most_constrained_negative, 0, core_member_a),
#'     core_member_b_negative = ifelse(
#'       share_trade_club_negative == 
#'         most_constrained_negative, 0, core_member_b)) %>% 
#'   mutate_at(vars(core_member_a_independent:core_member_b_negative),
#'             ~ifelse(is.na(.), 0, .)) %>% 
#'   # Now re-calculate the trade to club variables 
#'   mutate(trade_two_club_independent = ifelse(core_member_b == 1, 
#'                                              share_dyadic_independent, NA),
#'          trade_two_club_positive = ifelse(core_member_b == 1, 
#'                                           share_dyadic_positive, NA),
#'          trade_two_club_negative = ifelse(core_member_b == 1, 
#'                                           share_dyadic_negative, NA)) %>% 
#'   # Get share of trade to all club members
#'   group_by(id_a) %>% 
#'   mutate(share_trade_two_club_independent = sum(trade_two_club_independent, na.rm = T),
#'          share_trade_two_club_positive = sum(trade_two_club_positive, na.rm = T),
#'          share_trade_two_club_negative = sum(trade_two_club_negative, na.rm = T)) %>% 
#'   ungroup() %>% 
#'   # filter(id_a %in% c(1, 110, 120, 140),
#'   #        id_b %in% c(1, 110, 120, 140)) %>% 
#'   # Define lambda as second most constrained
#'   mutate(lambda_two_independent = min(second_constrained_independent, na.rm = T),
#'          lambda_two_positive = min(second_constrained_positive, na.rm = T),
#'          lambda_two_negative = min(second_constrained_negative, na.rm = T)) %>% 
#'   # Calculate change in lambda 
#'   mutate(delta_lambda_independent = second_constrained_independent - most_constrained_independent,
#'          delta_lambda_positive = second_constrained_positive - most_constrained_positive,
#'          delta_lambda_negative = second_constrained_negative - most_constrained_negative) %>% 
#'   distinct(climate_ideal_a, theta, share_trade_two_club_independent, 
#'            share_trade_two_club_positive, share_trade_two_club_negative,
#'            lambda_two_independent, lambda_two_positive, lambda_two_negative,
#'            delta_lambda_independent, delta_lambda_positive, delta_lambda_negative, ) %>% 
#'   # Calculate utilities to membership
#'   mutate(beta = 1) %>% # Define weights for climate versus trade
#'   rowwise() %>% 
#'   # Calculate utilities for the *independent* world
#'   # Utility to join
#'   mutate(utility_join_independent = (-abs(climate_ideal_a - theta)) + # Climate
#'            (1 - share_trade_two_club_independent) * (1 - lambda_two_independent)*beta + # Non-club trade
#'            (share_trade_two_club_independent*beta) - # Club trade
#'            (theta - lambda_two_independent)) %>%  # Aversion
#'   # Utility to abstain
#'   mutate(utility_abstain_independent = (-climate_ideal_a) + # Climate
#'            (1-share_trade_two_club_independent)*beta + # Non-club trade
#'            (share_trade_two_club_independent*(1-lambda_two_independent)*beta) + # Club trade
#'            (0)) %>%  # Aversion
#'   # Net, membership
#'   mutate(net_utility_independent = utility_join_independent - utility_abstain_independent,
#'          membership_independent = ifelse(net_utility_independent >= 0, 1, 0)) %>% 
#'   # Calculate utilities for the *positive* world
#'   # Utility to join
#'   mutate(utility_join_positive = (-abs(climate_ideal_a - theta)) + # Climate
#'            (1 - share_trade_two_club_positive) * (1 - lambda_two_positive)*beta + # Non-club trade
#'            (share_trade_two_club_positive*beta) - # Club trade
#'            (theta - lambda_two_positive)) %>%  # Aversion
#'   # Utility to abstain
#'   mutate(utility_abstain_positive = (-climate_ideal_a) + # Climate
#'            (1-share_trade_two_club_positive)*beta + # Non-club trade
#'            (share_trade_two_club_positive*(1-lambda_two_positive)*beta) + # Club trade
#'            (0)) %>%  # Aversion
#'   # Net, membership
#'   mutate(net_utility_positive = utility_join_positive - utility_abstain_positive,
#'          membership_positive = ifelse(net_utility_positive >= 0, 1, 0)) %>% 
#'   # Calculate utilities for the *negative* world
#'   # Utility to join
#'   mutate(utility_join_negative = (-abs(climate_ideal_a - theta)) + # Climate
#'            (1 - share_trade_two_club_negative) * (1 - lambda_two_negative)*beta + # Non-club trade
#'            (share_trade_two_club_negative*beta) - # Club trade
#'            (theta - lambda_two_negative)) %>%  # Aversion
#'   # Utility to abstain
#'   mutate(utility_abstain_negative = (-climate_ideal_a) + # Climate
#'            (1-share_trade_two_club_negative)*beta + # Non-club trade
#'            (share_trade_two_club_negative*(1-lambda_two_negative)*beta) + # Club trade
#'            (0)) %>%  # Aversion
#'   # Net, membership
#'   mutate(net_utility_negative = utility_join_negative - utility_abstain_negative,
#'          membership_negative = ifelse(net_utility_negative >= 0, 1, 0))